Spatial Transcriptomics Exploration

In this notebook we are comparing the use of Alevin-fry and Spaceranger for quantifying spatial transcriptomics libraries. Two spatial transcriptomic libraries were quantified using Alevin-fry and Spaceranger and the results were combined following the Alevin-fry tutorial. We are comparing that to if we were to not integrate the Spaceranger data with Alevin-fry and only use Spaceranger.

Additionally in this notebook we also explore the use of storing spatial data in a SpatialExperiment object or in a Seurat object.

Set Up

library(magrittr)
library(ggplot2)
library(SingleCellExperiment)
library(SpatialExperiment)
library(ggspavis)

Registered S3 method overwritten by 'ggside':
  method from   
  +.gg   ggplot2

library(ggupset)

# load in benchmarking functions that will be used for copying data and generating sample tables
function_path <- file.path(".." ,"benchmarking-functions", "R")
file.path(function_path, list.files(function_path, pattern = "*.R$")) %>%
  purrr::walk(source)

# set up file paths 
base_dir <- here::here()

# output folder to store alevin-fry and cellranger quants from S3
data_dir <- file.path(base_dir, "data", "spatial") 
s3_out <- file.path(data_dir, "data", "quants")

# output folder for image files 
image_dir <- file.path(data_dir, "data", "images")

# results directory 
results_dir <- file.path(data_dir, "results")

# create directories
if(!dir.exists(data_dir)){
  dir.create(data_dir, recursive = TRUE)
}

if(!dir.exists(image_dir)){
  dir.create(image_dir, recursive = TRUE)
}

if(!dir.exists(results_dir)){
  dir.create(results_dir, recursive = TRUE)
}

mito_file <- file.path(base_dir, "sample-info", "Homo_sapiens.GRCh38.103.mitogenes.txt")
  
# read in mito genes 
mito_genes <- readr::read_tsv(mito_file, col_names = "gene_id")

Rows: 111 Columns: 1
── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr (1): gene_id

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

mito_genes <- mito_genes %>%
  dplyr::pull(gene_id) %>%
  unique()

Copy data from S3

Comparing SpatialExperiment to Seurat

Importing Alevin-fry and Spaceranger output into a SpatialExperiment

First we are testing use of a SpatialExperiment object to store the spatial transcriptomics results. We will create two functions, one to deal with creating the combined object from both Alevin-fry and Spaceranger output, and the second to import Spaceranger output only.

# custom function for converting Alevin-fry and Spaceranger output to spe
convert_fry_spaceranger_spe <- function(fry_dir, spaceranger_dir, sample_name){
  fry_sce <- scpcaTools::read_alevin(fry_dir, 
                                     usa_mode = TRUE,
                                     which_counts = "spliced")
  
  # read in image data
  image_data <- SpatialExperiment::readImgData(spaceranger_dir,
                            sample_id = sample_name)
  
  # read in spatial coordinates 
  spatial_data_file <- file.path(spaceranger_dir, "tissue_positions_list.csv")
  spatial_data <- readr::read_csv(spatial_data_file, col_names = c("barcode", "in_tissue", "array_row",
                                                                   "array_col","pxl_row_in_fullres",
                                                                   "pxl_col_in_fullres"))
  # find common cells between Alevin-fry and spaceranger 
  spatial_data <- spatial_data %>%
    dplyr::mutate(barcode = gsub("-1", "", barcode))
  common_cells <- intersect(colnames(fry_sce), spatial_data$barcode)
  
  # subset fry_sce by common cells
  fry_sce_subset <- fry_sce[,common_cells]
  spatial_data_subset <- spatial_data %>%
    dplyr::filter(barcode %in% common_cells)
  
  # construct 'SpatialExperiment'
  fry_spe <- SpatialExperiment(
    assays = list(counts = assay(fry_sce_subset)),
    colData = colData(fry_sce_subset), 
    imgData = image_data,
    spatialData = DataFrame(spatial_data_subset),
    spatialCoordsNames = c("pxl_col_in_fullres", "pxl_row_in_fullres"),
    sample_id = sample_name)
}

spaceranger_only_spe <- function(spaceranger_outs_dir, sample_name){
  # read in spaceranger outputs directly using read10XVisium
  spaceranger_spe <- read10xVisium(spaceranger_outs_dir,
                                 sample_id = sample_name,
                                 type = "sparse",
                                 data = "raw",
                                 images = "lowres", 
                                 load = FALSE)
}

Let’s test these functions out and read in the Alevin-fry and Spaceranger output for one sample.

# get path to fry output directory 
fry_dir <- file.path(s3_out, "alevin-fry-knee", sample_ids[1])
fry_dir <- paste0(fry_dir, "-spliced_intron_txome_k31-salign-cr-like-em-knee")

# paths to spatial folders 
cellranger_folders <- paste0(sample_ids, "-cdna-spatial")
spatial_dir <- file.path(s3_out, "cellranger", cellranger_folders[1], "outs", "spatial")

# read in combined fry and spaceranger spe 
fry_spe <- convert_fry_spaceranger_spe(fry_dir, 
                                       spatial_dir, 
                                       sample_ids[1])

Rows: 4992 Columns: 6
── Column specification ──────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): barcode
dbl (5): in_tissue, array_row, array_col, pxl_row_in_fullres, pxl_col_in_fullres

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’

Below we import just the Spaceranger outputs into a SpatialExperiment object.

# spaceranger output paths
spaceranger_out <- file.path(s3_out, "cellranger", cellranger_folders[1], "outs")

# read in spaceranger output directly using read10XVisium
spaceranger_spe <- spaceranger_only_spe(spaceranger_out, sample_ids[1])

Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’

Importing as Seurat Object

Below we will import the same data but as a Seurat object instead of a SpatialExperiment object.

# first create sce from Alevin output
fry_sce <- scpcaTools::read_alevin(fry_dir,
                                   usa_mode = TRUE,
                                   which_counts = "spliced")

# convert sce to seurat object 
fry_seurat <- CreateSeuratObject(counts = counts(fry_sce),
                                 project = "SPATIAL",
                                 assay = "Spatial")

# read in image data using Seurat and subset to common cells
image_data <- Read10X_Image(spatial_dir)
rownames(image_data@coordinates) <- gsub("-1","", rownames(image_data@coordinates))
common_cells <- intersect(Cells(x = fry_seurat), rownames(image.data@coordinates))

Error in h(simpleError(msg, call)) : 
  error in evaluating the argument 'y' in selecting a method for function 'intersect': error in evaluating the argument 'x' in selecting a method for function 'rownames': object 'image.data' not found

# load spaceranger outputs directly into Seurat object 
spaceranger_seurat <- Load10X_Spatial(spaceranger_out)

'giveCsparse' has been deprecated; setting 'repr = "T"' for you

Comparing Seurat and Spatial Experiment

Now let’s take a look at how each of these objects stores information for the experiment. We will just look at the files that have been produced from Spaceranger only and compare the Seurat object to the SpatialExperiment.

head(colData(spaceranger_spe))

DataFrame with 6 rows and 1 column
                     sample_id
                   <character>
AAACAACGAATAGTTC-1 SCPCR000372
AAACAAGTATCTCCCA-1 SCPCR000372
AAACAATCTACTAGCA-1 SCPCR000372
AAACACCAATAACTGC-1 SCPCR000372
AAACAGAGCGACTCCT-1 SCPCR000372
AAACAGCTTTCAGAAG-1 SCPCR000372

head(spaceranger_seurat@meta.data)

Seurat automatically adds some per cell QC metrics here, but in general the colData looks similar to what we would expect for a SingleCellExperiment for the SpatialExperiment.

head(rowData(spaceranger_spe))

DataFrame with 6 rows and 1 column
                     symbol
                <character>
ENSG00000186092       OR4F5
ENSG00000284733      OR4F29
ENSG00000284662      OR4F16
ENSG00000187634      SAMD11
ENSG00000188976       NOC2L
ENSG00000187961      KLHL17

head(spaceranger_seurat[["Spatial"]]@meta.features)

spaceranger_spe <- spaceranger_spe %>%
  scuttle::addPerCellQCMetrics(subsets = list(mito_genes = mito_genes[mito_genes %in% rownames(spaceranger_spe)])) 

Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’

head(colData(spaceranger_spe))

DataFrame with 6 rows and 7 columns
                     sample_id       sum  detected subsets_mito_genes_sum subsets_mito_genes_detected
                   <character> <numeric> <integer>              <numeric>                   <integer>
AAACAACGAATAGTTC-1 SCPCR000372      9094      3963                    148                          10
AAACAAGTATCTCCCA-1 SCPCR000372      7070      2993                    304                          10
AAACAATCTACTAGCA-1 SCPCR000372     14917      5212                    313                          11
AAACACCAATAACTGC-1 SCPCR000372      3999      2280                     93                          10
AAACAGAGCGACTCCT-1 SCPCR000372     14765      5050                    139                          10
AAACAGCTTTCAGAAG-1 SCPCR000372      3972      2174                    134                           9
                   subsets_mito_genes_percent     total
                                    <numeric> <numeric>
AAACAACGAATAGTTC-1                   1.627447      9094
AAACAAGTATCTCCCA-1                   4.299859      7070
AAACAATCTACTAGCA-1                   2.098277     14917
AAACACCAATAACTGC-1                   2.325581      3999
AAACAGAGCGACTCCT-1                   0.941416     14765
AAACAGCTTTCAGAAG-1                   3.373615      3972

You can still use the same functions to add in per cell metrics that we are using in scpcaTools now on a SpatialExperiment.

Let’s specifically look at how each of them stores the data related to the spatial information. The SpatialExperiment has three places where it stores spatial information, spatialData, spatialCoords, and imgData.

spatialData(spaceranger_spe)

DataFrame with 4992 rows and 3 columns
                   in_tissue array_row array_col
                   <logical> <integer> <integer>
AAACAACGAATAGTTC-1     FALSE         0        16
AAACAAGTATCTCCCA-1      TRUE        50       102
AAACAATCTACTAGCA-1     FALSE         3        43
AAACACCAATAACTGC-1     FALSE        59        19
AAACAGAGCGACTCCT-1      TRUE        14        94
...                      ...       ...       ...
TTGTTTCACATCCAGG-1      TRUE        58        42
TTGTTTCATTAGTCTA-1      TRUE        60        30
TTGTTTCCATACAACT-1      TRUE        45        27
TTGTTTGTATTACACG-1      TRUE        73        41
TTGTTTGTGTAAATTC-1      TRUE         7        51

head(spatialCoords(spaceranger_spe))

                   pxl_col_in_fullres pxl_row_in_fullres
AAACAACGAATAGTTC-1               1370               2280
AAACAAGTATCTCCCA-1               6047               6970
AAACAATCTACTAGCA-1               1642               3744
AAACACCAATAACTGC-1               6927               2483
AAACAGAGCGACTCCT-1               2659               6513
AAACAGCTTTCAGAAG-1               5424               1930

spatialData holds information about whether or not each barcode was covered by tissue and spatial metadata. spatialCoords holds the coordinates for each spot.

This same information can be found by looking in the images slot in the Seurat object.

imgData holds information about each image that is stored as multiple slices of an image can be stored within one SpatialExperiment.

imgData(spaceranger_spe)

DataFrame with 1 row and 4 columns
    sample_id    image_id   data scaleFactor
  <character> <character> <list>   <numeric>
1 SCPCR000372      lowres   ####        0.06

This same information can be found by looking at the images slot in the Seurat object.

spaceranger_seurat@images$slice1@scale.factors

$spot
[1] 0.2

$fiducial
[1] 113.7005

$hires
[1] 0.2

$lowres
[1] 0.06

attr(,"class")
[1] "scalefactors"

Let’s take a look at some example plots that can be produced with each of these objects.

# plot tissue spots coloring by total UMI count
plotSpots(spe = spaceranger_spe, 
          x_coord = "pxl_col_in_fullres", 
          y_coord = "pxl_row_in_fullres", 
          annotate = "sum") +
  scale_color_viridis_c()

Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing
scale.

SpatialFeaturePlot(spaceranger_seurat, features = "nCount_Spatial")

`guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.

It looks like overall they both would work and that the SpatialExperiment will provide similar functionality to how we have things already set up for single-cell experiments. I think I have a slight preference for the Seurat plots as you can adjust the transparency of the points and also see the tissue, but I also think we may be able to play around with the plots with the SpatialExperiment more if we want something similar to that.

Alevin-fry + Spaceranger versus Spaceranger Only

Now let’s take a look at comparing the two methods of using Alevin-fry + Spaceranger to only Spaceranger for quantification. To do this, we will read in a list of Alevin-fry + Spaceranger SpatialExperiment objects and Spaceranger objects and then merge them into one list before grabbing the per cell and per gene quality metrics.

Create Spatial Experiments

# create sample info dataframe to be joined with per cell dataframe later
sample_info_df <- quant_info_table(data_dir= s3_out, 
                 tools = c("cellranger", "alevin-fry-knee"),
                 samples = sample_ids) %>%
  # convert cellranger to spaceranger 
  dplyr::mutate(tool = ifelse(tool == "cellranger", "spaceranger", tool))

# modify sample data frame to get spatial and fry directories for each sample
input_df <- sample_info_df %>%
  dplyr::select(tool, quant_dir, sample, data_dir) %>%
  dplyr::mutate(data_dir = dplyr::case_when(tool == "spaceranger" ~ 
                                            file.path(data_dir, "outs", "spatial"),
                                            tool == "alevin-fry" ~ data_dir)) %>%
  tidyr::pivot_wider(id_cols = sample, names_from = tool, values_from = data_dir)

# use custome function to create list of spe's
fry_spe_list <- mapply(convert_fry_spaceranger_spe,
                       input_df$`alevin-fry`,
                       input_df$spaceranger,
                       input_df$sample)

Rows: 4992 Columns: 6
── Column specification ──────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): barcode
dbl (5): in_tissue, array_row, array_col, pxl_row_in_fullres, pxl_col_in_fullres

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’
Rows: 4992 Columns: 6
── Column specification ──────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): barcode
dbl (5): in_tissue, array_row, array_col, pxl_row_in_fullres, pxl_col_in_fullres

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’

names(fry_spe_list) <- paste0(sample_ids, "-fry-spaceranger-combined")

# create list of spaceranger only processed spe's 
# spaceranger output paths
spaceranger_outs <- file.path(s3_out, "cellranger", cellranger_folders, "outs")
names(spaceranger_outs) <- sample_ids

spaceranger_spe_list <- purrr::imap(spaceranger_outs, spaceranger_only_spe)

Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’

names(spaceranger_spe_list) <- paste0(sample_ids, "-spaceranger-only")

all_spe_list <- append(fry_spe_list, spaceranger_spe_list)

# calculate per cell QC and output to a combined data frame with plotting 
all_spe_list <- all_spe_list %>%
    purrr::map(
      ~ scuttle::addPerCellQCMetrics(.x, 
                                     subsets = list(mito = mito_genes[mito_genes %in% rownames(.x)]))) %>%
  purrr::map(scuttle::addPerFeatureQCMetrics)

Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’

Per Spot QC Metrics

# create custom function to grab colData and convert to data frame 
spatial_coldata_to_df <- function(spe){
  df <- as.data.frame(colData(spe)) %>%
    tibble::rownames_to_column(var = "spot_id")
}

fry_samples = c("SCPCR000372-fry-spaceranger-combined", "SCPCR000373-fry-spaceranger-combined")
spaceranger_samples = c("SCPCR000372-spaceranger-only", "SCPCR000373-spaceranger-only")

# join coldata dataframe with sample info
coldata_df <- purrr::map_df(all_spe_list,spatial_coldata_to_df, .id = "sample") %>%
  dplyr::mutate(tool = dplyr::case_when(sample %in% fry_samples ~ "alevin-fry",
                                        sample %in% spaceranger_samples ~ "spaceranger"),
                # extract sample name only to allow for joining
                sample = stringr::word(sample, 1, sep = "-"),
                # remove extra -1 from spaceranger barcodes
                spot_id = gsub("-1", "", spot_id))%>%
  dplyr::left_join(sample_info_df,
                   by = c("sample", "tool"))

# identify shared spots only 
spot_counts <- coldata_df %>%  
  dplyr::count(spot_id, sample)

common_spots <- spot_counts %>%
  dplyr::filter(n == 2) %>%
  dplyr::pull(spot_id)

coldata_df_common <- coldata_df %>%
  dplyr::filter(spot_id %in% common_spots)

In addition to filtering on shared spots we will also want to filter based on spots that are found to overlap with tissue.

# also want to filter the spe's directly 
filter_spe <- function(spe){
  spe <- spe[, spatialData(spe)$in_tissue == 1]
}

all_spes_filter <- all_spe_list %>%
  purrr::map(filter_spe)

Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’

# custom function for plotting spe results 
plot_spe <- function(spe, sample, column){
  #column_sym <- rlang::sym(column)
  plotSpots(spe, 
            x_coord = "pxl_col_in_fullres", 
            y_coord = "pxl_row_in_fullres", 
            annotate = column) +
  scale_color_viridis_c() + 
    ggtitle(sample)
}

First we will look at the per cell metrics: mitochondrial reads per cell, total UMI per cell, and total genes detected per cell.

# % mitochondrial reads/ spot 
ggplot(coldata_df_common, aes(x = tool, y = subsets_mito_percent, fill = tool)) + 
  geom_boxplot() + 
  facet_grid(~ sample) +
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  ylab("Mito Percent") + 
  xlab("")

all_spes_filter %>%
  purrr::imap(plot_spe, "subsets_mito_percent")

Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing
scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing
scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing
scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing
scale.

$`SCPCR000372-fry-spaceranger-combined`

$`SCPCR000373-fry-spaceranger-combined`

$`SCPCR000372-spaceranger-only`

$`SCPCR000373-spaceranger-only`

Overall it looks like mitochondrial content is low and fairly similar across both tools.

all_spes_filter %>%
  purrr::imap(plot_spe, "sum")

Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing
scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing
scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing
scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing
scale.

$`SCPCR000372-fry-spaceranger-combined`

$`SCPCR000373-fry-spaceranger-combined`

$`SCPCR000372-spaceranger-only`

$`SCPCR000373-spaceranger-only`

Interestingly when looking at the pictures (besides that the Alevin-fry and Spaceranger pictures are flipped), there appears to be different patterns of total UMI count per spot showing up between fry and Spaceranger.

all_spes_filter %>%
  purrr::imap(plot_spe, "detected")

Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing
scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing
scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing
scale.
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing
scale.

$`SCPCR000372-fry-spaceranger-combined`

$`SCPCR000373-fry-spaceranger-combined`

$`SCPCR000372-spaceranger-only`

$`SCPCR000373-spaceranger-only`

Generally it looks like both tools are fairly similar in the number of UMI/cell and genes detected/cell but the actual spatial distribution of looks very different within each sample across both tools which is concerning.

Per Gene QC Metrics

Let’s also look at the correlation of mean gene expression across shared genes.

spatial_rowdata_to_df <- function(spe){
  df <- as.data.frame(rowData(spe)) %>%
    tibble::rownames_to_column(var = "gene_id")
}

# combine rowdata with sample info
rowdata_df <- purrr::map_df(all_spe_list,spatial_rowdata_to_df, .id = "sample") %>%
  dplyr::mutate(tool = dplyr::case_when(sample %in% fry_samples ~ "alevin-fry",
                                        sample %in% spaceranger_samples ~ "spaceranger"),
                sample = stringr::word(sample, 1, sep = "-")) %>%
  dplyr::left_join(sample_info_df,
                   by = c("sample", "tool"))


gene_counts <- rowdata_df %>% 
  # remove genes that have a low frequency of being detected
  dplyr::filter(detected >= 5.0) %>%
  dplyr::count(gene_id, sample)

# restrict to only common genes 
common_genes <- gene_counts %>%
  dplyr::filter(n == 2) %>%
  dplyr::pull(gene_id)

rowdata_df_common <- rowdata_df %>%
  dplyr::filter(
    (gene_id %in% common_genes) 
  )

rowdata_cor <- rowdata_df_common %>%
  dplyr::select(tool, gene_id, sample, mean) %>%
  # spread the mean expression stats to one column per caller
  tidyr::pivot_wider(id_cols = c(gene_id, sample),
                     names_from = c("tool"),
                     values_from = mean) %>%
  # drop rows with NA values to ease correlation calculation
  tidyr::drop_na()

# look at correlation between the two tools
rowdata_cor %>% 
  dplyr::group_by(sample) %>%
  dplyr::summarize(
    alevin_fry_knee_spaceranger_cor = cor(`spaceranger`, `alevin-fry`, method = "spearman")
  )

# mean gene expression across shared genes 
ggplot(rowdata_cor, aes(x = `spaceranger`, y = `alevin-fry`)) +
  geom_point(size = 0.5, alpha = 0.1) + 
  facet_wrap(~ sample) + 
  scale_x_log10() + 
  scale_y_log10() + 
  labs(x = "Spaceranger mean gene expression", y = "Alevin Fry Mean gene expression") + 
  theme_classic()

Shared genes across tools

The final thing we will look at is the overlap of genes detected in each tool. First we need to filter the SpatialExperiment objects to only have spots that are found in both tools and then we filter to include only genes that are found in both indices before looking at the overlap.

# make a function to filter sces by subset of spots and re-calculate feature stats
filter_spe <- function(spe, spots){
  # remove "-1" at end of barcode for spaceranger spes
  colnames(spe) <- gsub("-1", "", colnames(spe))
  cells_to_keep <- colnames(spe) %in% spots
  rowData(spe) <- NULL
  spe[, cells_to_keep] %>%
    scuttle::addPerFeatureQCMetrics()
}

# filter all spes to only have common spots 
spe_list_common <- all_spes_filter %>%
  purrr::map(filter_spe, spots = common_spots)

Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’
Found more than one class "SpatialImage" in cache; using the first, from namespace 'SeuratObject'
Also defined by ‘SpatialExperiment’

# grab rowdata from filtered sces 
rowdata_df_filtered <- purrr::map_df(spe_list_common, spatial_rowdata_to_df, .id = "sample") %>%
  dplyr::mutate(tool = dplyr::case_when(sample %in% fry_samples ~ "alevin-fry",
                                        sample %in% spaceranger_samples ~ "spaceranger"),
                sample = stringr::word(sample, 1, sep = "-")) %>%
  dplyr::left_join(sample_info_df,
                   by = c("sample", "tool"))

# get genes common in all tools
common_genes <- rowdata_df_filtered %>%
  dplyr::select(gene_id, tool) %>%
  dplyr::distinct() %>%
  dplyr::group_by(gene_id) %>%
  dplyr::tally() %>%
  # filter for genes found in both spaceranger and alevin-fry
  dplyr::filter(n == 2) %>%
  dplyr::pull(gene_id)

# filter rowdata_df to only include genes found in all tools and genes with mean > 0 and detected > 0 in all cells
rowdata_df_filtered <- rowdata_df_filtered %>% 
  dplyr::filter(gene_id %in% common_genes) %>%
  dplyr::filter(mean > 0 & detected > 0)

# create upset plot 
gene_detect_df <- rowdata_df_filtered %>%
  dplyr::group_by(gene_id) %>%
  dplyr::summarise(tools_detected = list(unique(tool)))

ggplot(gene_detect_df, aes(x = tools_detected)) +
  geom_bar() +
  scale_x_upset(n_intersections = 3)

Similar to with the single-cell data most genes are detected in both Alevin-fry and Spaceranger which is good.

Concluding thoughts

• It appears that we can use the SpatialExperiment to store the quantification output and it will maintain similar functionality to the SingleCellExperiment obejcts we use for single-cell and single-nuclei data.
• Although Alevin-fry and Spaceranger only result in similar distributions of UMI/spot and genes/spot, the patterning across the tissue is not the same which is not ideal and doesn’t provide support for adding the use of Alevin-fry to our workflow (since we have to use Spaceranger anyways).

Session Info

sessioninfo::session_info()

─ Session info  🥦  🧔🏿  😋   ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
 hash: broccoli, person: dark skin tone, beard, face savoring food

 setting  value
 version  R version 4.1.1 (2021-08-10)
 os       macOS Catalina 10.15.7
 system   x86_64, darwin17.0
 ui       RStudio
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       America/Chicago
 date     2021-11-17
 rstudio  1.4.1106 Tiger Daylily (desktop)
 pandoc   2.11.4 @ /Applications/RStudio.app/Contents/MacOS/pandoc/ (via rmarkdown)

─ Packages ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
 package              * version  date (UTC) lib source
 assertthat             0.2.1    2019-03-21 [1] CRAN (R 4.1.0)
 beachmat               2.10.0   2021-10-26 [1] Bioconductor
 Biobase              * 2.54.0   2021-10-26 [1] Bioconductor
 BiocGenerics         * 0.40.0   2021-10-26 [1] Bioconductor
 BiocParallel           1.28.0   2021-10-26 [1] Bioconductor
 bit                    4.0.4    2020-08-04 [1] CRAN (R 4.1.0)
 bit64                  4.0.5    2020-08-30 [1] CRAN (R 4.1.0)
 bitops                 1.0-7    2021-04-24 [1] CRAN (R 4.1.0)
 bslib                  0.3.1    2021-10-06 [1] CRAN (R 4.1.0)
 cli                    3.1.0    2021-10-27 [1] CRAN (R 4.1.1)
 colorspace             2.0-2    2021-06-24 [1] CRAN (R 4.1.0)
 crayon                 1.4.2    2021-10-29 [1] CRAN (R 4.1.0)
 data.table             1.14.2   2021-09-27 [1] CRAN (R 4.1.0)
 DBI                    1.1.1    2021-01-15 [1] CRAN (R 4.1.0)
 DelayedArray           0.20.0   2021-10-26 [1] Bioconductor
 DelayedMatrixStats     1.16.0   2021-10-26 [1] Bioconductor
 digest                 0.6.28   2021-09-23 [1] CRAN (R 4.1.0)
 dplyr                  1.0.7    2021-06-18 [1] CRAN (R 4.1.0)
 dqrng                  0.3.0    2021-05-01 [1] CRAN (R 4.1.0)
 DropletUtils           1.14.1   2021-11-08 [1] Bioconductor
 edgeR                  3.36.0   2021-10-26 [1] Bioconductor
 ellipsis               0.3.2    2021-04-29 [1] CRAN (R 4.1.0)
 evaluate               0.14     2019-05-28 [1] CRAN (R 4.1.0)
 fansi                  0.5.0    2021-05-25 [1] CRAN (R 4.1.0)
 fastmap                1.1.0    2021-01-25 [1] CRAN (R 4.1.0)
 generics               0.1.1    2021-10-25 [1] CRAN (R 4.1.0)
 GenomeInfoDb         * 1.30.0   2021-10-26 [1] Bioconductor
 GenomeInfoDbData       1.2.7    2021-11-16 [1] Bioconductor
 GenomicRanges        * 1.46.0   2021-10-26 [1] Bioconductor
 ggplot2              * 3.3.5    2021-06-25 [1] CRAN (R 4.1.0)
 ggside                 0.1.3    2021-10-24 [1] CRAN (R 4.1.0)
 ggspavis             * 1.0.0    2021-10-26 [1] Bioconductor
 ggupset              * 0.3.0    2020-05-05 [1] CRAN (R 4.1.0)
 glue                   1.5.0    2021-11-07 [1] CRAN (R 4.1.0)
 gtable                 0.3.0    2019-03-25 [1] CRAN (R 4.1.0)
 HDF5Array              1.22.1   2021-11-14 [1] Bioconductor
 here                   1.0.1    2020-12-13 [1] CRAN (R 4.1.0)
 hms                    1.1.1    2021-09-26 [1] CRAN (R 4.1.0)
 htmltools              0.5.2    2021-08-25 [1] CRAN (R 4.1.0)
 htmlwidgets            1.5.4    2021-09-08 [1] CRAN (R 4.1.0)
 IRanges              * 2.28.0   2021-10-26 [1] Bioconductor
 jquerylib              0.1.4    2021-04-26 [1] CRAN (R 4.1.0)
 jsonlite               1.7.2    2020-12-09 [1] CRAN (R 4.1.0)
 knitr                  1.36     2021-09-29 [1] CRAN (R 4.1.0)
 lattice                0.20-45  2021-09-22 [1] CRAN (R 4.1.0)
 lifecycle              1.0.1    2021-09-24 [1] CRAN (R 4.1.0)
 limma                  3.50.0   2021-10-26 [1] Bioconductor
 locfit                 1.5-9.4  2020-03-25 [1] CRAN (R 4.1.0)
 magick                 2.7.3    2021-08-18 [1] CRAN (R 4.1.0)
 magrittr             * 2.0.1    2020-11-17 [1] CRAN (R 4.1.0)
 Matrix                 1.3-4    2021-06-01 [1] CRAN (R 4.1.1)
 MatrixGenerics       * 1.6.0    2021-10-26 [1] Bioconductor
 matrixStats          * 0.61.0   2021-09-17 [1] CRAN (R 4.1.0)
 munsell                0.5.0    2018-06-12 [1] CRAN (R 4.1.0)
 pillar                 1.6.4    2021-10-18 [1] CRAN (R 4.1.1)
 pkgconfig              2.0.3    2019-09-22 [1] CRAN (R 4.1.0)
 purrr                  0.3.4    2020-04-17 [1] CRAN (R 4.1.0)
 R.methodsS3            1.8.1    2020-08-26 [1] CRAN (R 4.1.0)
 R.oo                   1.24.0   2020-08-26 [1] CRAN (R 4.1.0)
 R.utils                2.11.0   2021-09-26 [1] CRAN (R 4.1.0)
 R6                     2.5.1    2021-08-19 [1] CRAN (R 4.1.0)
 Rcpp                   1.0.7    2021-07-07 [1] CRAN (R 4.1.0)
 RCurl                  1.98-1.5 2021-09-17 [1] CRAN (R 4.1.0)
 readr                  2.1.0    2021-11-11 [1] CRAN (R 4.1.0)
 rhdf5                  2.38.0   2021-10-26 [1] Bioconductor
 rhdf5filters           1.6.0    2021-10-26 [1] Bioconductor
 Rhdf5lib               1.16.0   2021-10-26 [1] Bioconductor
 rjson                  0.2.20   2018-06-08 [1] CRAN (R 4.1.0)
 rlang                  0.4.12   2021-10-18 [1] CRAN (R 4.1.1)
 rmarkdown              2.11     2021-09-14 [1] CRAN (R 4.1.0)
 rprojroot              2.0.2    2020-11-15 [1] CRAN (R 4.1.0)
 rstudioapi             0.13     2020-11-12 [1] CRAN (R 4.1.0)
 S4Vectors            * 0.32.2   2021-11-07 [1] Bioconductor
 sass                   0.4.0    2021-05-12 [1] CRAN (R 4.1.0)
 scales                 1.1.1    2020-05-11 [1] CRAN (R 4.1.0)
 scpcaTools             0.1.2    2021-10-12 [1] Github (AlexsLemonade/scpcaTools@2cdad4c)
 scuttle                1.4.0    2021-10-26 [1] Bioconductor
 sessioninfo            1.2.1    2021-11-02 [1] CRAN (R 4.1.0)
 SeuratObject         * 4.0.3    2021-11-10 [1] CRAN (R 4.1.0)
 SingleCellExperiment * 1.16.0   2021-10-26 [1] Bioconductor
 sparseMatrixStats      1.6.0    2021-10-26 [1] Bioconductor
 SpatialExperiment    * 1.4.0    2021-10-26 [1] Bioconductor
 stringi                1.7.5    2021-10-04 [1] CRAN (R 4.1.1)
 stringr                1.4.0    2019-02-10 [1] CRAN (R 4.1.0)
 SummarizedExperiment * 1.24.0   2021-10-26 [1] Bioconductor
 tibble                 3.1.6    2021-11-07 [1] CRAN (R 4.1.0)
 tidyr                  1.1.4    2021-09-27 [1] CRAN (R 4.1.0)
 tidyselect             1.1.1    2021-04-30 [1] CRAN (R 4.1.0)
 tzdb                   0.2.0    2021-10-27 [1] CRAN (R 4.1.1)
 utf8                   1.2.2    2021-07-24 [1] CRAN (R 4.1.0)
 vctrs                  0.3.8    2021-04-29 [1] CRAN (R 4.1.0)
 vroom                  1.5.6    2021-11-10 [1] CRAN (R 4.1.0)
 withr                  2.4.2    2021-04-18 [1] CRAN (R 4.1.0)
 xfun                   0.28     2021-11-04 [1] CRAN (R 4.1.0)
 XVector                0.34.0   2021-10-26 [1] Bioconductor
 yaml                   2.2.1    2020-02-01 [1] CRAN (R 4.1.0)
 zlibbioc               1.40.0   2021-10-26 [1] Bioconductor

 [1] /Library/Frameworks/R.framework/Versions/4.1/Resources/library

────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────